Modelling and optimization of photon pair sources based on spontaneous parametric 

down-conversion 
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We address the problem of efficient modelling of photon pairs generated in spontaneous parametric 
down-conversion and coupled into single-mode fibers. It is shown that when the range of relevant 
transverse wave vectors is restricted by the pump and fiber modes, the computational complexity 
can be reduced substantially with the help of the paraxial approximation, while retaining the full 
spectral characteristics of the source. This approach can serve as a basis for efficient numerical 
calculations, or can be combined with analytically tractable approximations of the phase matching 
function. We introduce here a cosine-gaussian approximation of the phase matching function which 
works for a broader range of parameters than the gaussian model used previously. The developed 
modelling tools are used to evaluate characteristics of the photon pair sources such as the pair 
production rate and the spectral purity quantifying frequency correlations. Strategies to generate 
spectrally uncorrelated photons, necessary in multiphoton interference experiments, are analyzed 
with respect to trade-offs between parameters of the source. 

PACS numbers: 42.65.Lm,42.50.Dv,03.67.Bg 



I. INTRODUCTION 

Spontaneous parametric down conversion (SPDC) is a 
nonlinear process in which a pump photon interacting 
with a crystal decays into two daughter photons. The 
process has been successfully employed to demonstrate 
fundamental aspects of quantum mechanics such as the 
violation of Bell's inequalities [H, 0], and utilized in im- 
plementations of quantum teleportation [!, 0, [H , quan- 
tum cryptography jf|, linear optical quantum informa- 
tion processing p}, and other quantum-enhanced tech- 
nologies. 

Typically, photon pairs emerging from non-linear me- 
dia are described by a complicated spatio-temporal wave 
function that exhibits correlations in multiple degrees 
of freedom. In contrast, many applications of photon 
pairs require their preparation in single isolated spatio- 
temporal modes. Single spatial modes can be selected by 
coupling photons into single-mode fibers (SMFs) that in 
effect filter heavily the SPDC light. Furthermore, many 
protocols rely on interference between photons originat- 
ing from independent sources @, [§]. Although spatial 
modes are well defined by SMFs, the interference visi- 
bility may be compromised by undesirable spectral cor- 
relations within individual pairs. One way to tailor the 
spectral degree of freedom is to use narrowband interfer- 
ence filters, which is easy to implement in an experiment, 
but in consequence reduces the useful photon flux. An 
alternative approach is to adjust the setup parameters 
to enforce the source to produce spectrally uncorrelated 
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pairs (lflUJllialia. 

These issues bring the question of optimizing the use- 
ful fraction of photon pairs produced by SPDC sources. 
A purely experimental approach would be just to try var- 
ious alignments of the source. In practice, this strategy 
would be rather burdensome owing to the large num- 
ber of controllable parameters of the setup, their time- 
consuming adjustments, and long data acquisition times. 
A natural alternative is to resort to numerical modeling. 
This however presents its own challenges, as including all 
relevant degrees of freedom is computationally demand- 
ing. 

In this paper we discuss approximate methods that al- 
leviate the numerical load necessary to model faithfully 
realistic SPDC sources. Our approach is based on an 
observation that optical fibers collecting photons define 
a relatively narrow range of wave vectors that need to 
be included in calculations. This justifies applying the 
paraxial approximation, which makes a substantial por- 
tion of the problem tractable analytically. The paraxial 
approximation can be also combined with a simplification 
of the two-photon wave function to an analytically man- 
ageable form leading to closed formulas. We exploit these 
strategies to analyze the performance of SPDC sources 
in quantum information applications. 

Coupling of down-converted photons into SMFs has 
been a subject of a number of works, especially in the 
case of cw pumping. Kurtsiefer et al. gave a simple 
argument showing that careful matching of the SPDC 
output with the fiber modes increases the collection ef- 
ficiency. Mathematical models for a cw-pumped source 
has been derived and compared with experimental data 
in Refs. [H, EE, The collinear case has been ana- 
lyzed theoretically in Ref. [l^]. Dragan [T0| used a gaus- 
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sian approximation to model fiber-coupled sources. The 
counterintuitive scaling of the production rates with the 
crystal length has been pointed out by Lee et at 
and a detailed analysis of quasi-phase matched structures 
has been presented by Ljunggren and Tengner (2p| . l2lj . 
More recently, Ling et al. [22] provided a method to esti- 
mate the absolute emission rates for cw pumping. In the 
present paper, we concentrate on pulse-pumped SPDC 
sources and optimization of their performance parame- 
ters. 

Our numerical calculations incorporate the exact form 
of dispersion relations for the nonlinear medium and use 
a second-order expansion of phase mismatch in the trans- 
verse wave vectors of the SPDC photons. The modeling 
is based on two strategies. The first approach resorts 
to numerical means, but with minimized computational 
effort that will nevertheless deliver highly accurate re- 
sults in a broad range of parameters. This method has 
been used in Refs. [23, [24j to compare experimentally 
measured characteristics of down-conversion sources with 
theoretical predictions. The second approach will provide 
expressions for the biphoton wave function in a closed 
analytical form through a further approximation to the 
phase matching functions. This approach, which we will 
call the co sine- gaus sian approximation (CGA) is based 
on a more accurate analytically integrable model of the 
phase matching function than the gaussian model studied 
previously [Tol.Tlll. |25| . We compare both the approaches 
with direct numerical calculations when no paraxial ap- 
proximation is applied and all integrals are evaluated by 
numerical means. As an application of the developed 
tools, we discuss generation of spectrally uncorrelated 
photons in a type-I /J-barium borate (BBO) crystal. We 
consider here two strategies to reduce spectral correla- 
tion: one method is to adjust carefully the pump pulse 
and collection modes, while the other one is to restrict 
the spectrum of the generated photons with the help of 
interference filters. We compare source brightness that 
can be achieved using both methods and relate these re- 
sults to previous discussions [Io| . 

The paper is organized as follows. In Sec.[ll]we present 
the setup under consideration and derive the biphoton 
wave function in free space. Section [TTT] presents basic 
assumptions about propagation of a pump beam and out- 
put photons and the impact of spatial filtering imposed 
by SMFs. The biphoton wave function within the parax- 
ial approximation is derived. Next in Sec.[lV]we present 
the cosine gaussian approximation and apply it to derive 
an analytical formula for the wave function of a photon 
pair coupled into SMFs. The figures of merit are defined 
in Sec.|Vl and the approximation of perfect phase match- 
ing is used to gain some basic intuitions. Next in Sec. lVIl 
we compare the computational effort and applicability 
of developed methods. Finally, in Sec. IVIII we analyze 
strategies to reduce spectral correlations within photon 
pairs. 



II. TWO-PHOTON WAVE FUNCTION 

In the non-degenerate down-conversion process, the 
pump field, described by the positive-frequency part of 
the electric field E p + \r, t), interacts with quantized sig- 
nal and idler fields, whose creation-operator parts will 



be denoted as M \r,t) and E\ '(r,t). The interaction 
hamiltonian has the form of an integral over the volume 
V of the crystal [Ifj]: 
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= f°|^ J d 3 rE^(r,t)Ei-\r,t)Et\T,t) 
v 



+ H.c, (II.l) 

where eo is the vacuum permittivity and x' 2 ' denotes the 
second-order nonlinear susceptibility coefficient, approx- 
imated by a constant. We will assume that the nonlinear 
interaction is weak enough to neglect pump depletion and 
to justify the first order perturbation theory. We will fo- 
cus here on type-I phase-matching, when both the down- 
converted photons have the same polarization direction, 
perpendicular to that of the pump pulse. The case of 
type-II phase matching can be analyzed analogously. 

We will take the nonlinear crystal to be a thin slab of 
thickness L oriented perpendicular to z-axis and extend- 
ing from z = —L/2 to z = L/2, as illustrated in Fig. [TJ 
The pump pulse propagates along z-direction outside the 
crystal. Following Rubin et al. [23] we parameterize the 
waves using the frequencies co and the wave vector com- 
ponents ki. perpendicular to z. These quantities are pre- 
served at the crystal-free space interface. In order to cal- 
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FIG. 1: (Color online) The geometry of a photon pair source. 
A crystal exhibiting nonlinearlity is pumped by gaussian 
pulse. The generated light emerging at angles a s and ati is 
coupled into single mode optical fibers. 

ibrate the pump power, it will be convenient to introduce 
a normalized pump pulse amplitude A p (k p ±, uj p ) satisfy- 
ing J d 2 k p xdwp|^4p(kp_L, uj p )\ 2 = 1. We will assume the 
pump pulse amplitude in a factorable form, with no spa- 
tiotemporal correlations: 



A P em P(iv)A s e(k ± ) 



(H.2) 
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where A p emp (w) refers to temporal and A p p (k±) to spatial 
part. Both parts will be taken in a gaussian form: 



A* emp (a 



■ cxp 



(u - 2w f 



A p *(k ± ) = -^expl -^ki 



(H.3) 



(II.4) 



where r p stands for the pulse duration, w p for the pump 
beam width, and 2wo is the central frequency of the pump 
pulse. 

The positive-frequency part of the pump pulse electric 
field Ep + \r,t) is the Fourier transform of the spectral 
amplitude: 

E ( p +\v,t)=£ p J d 2 Vdcj p A p (V,w p )e i ( k ^ r — *> 

(II.5) 

where £ p characterizes the strength of the pump pulse 
and the squared modulus \£ p \ 2 proportional to the pump 
pulse energy. Subsequently, we assume the following 
modal expansion for the signal s and idler i field op- 
erators: 

£<-)(r, t)=£ p J d 2 k M ±<H e- ik " r+i ^ t at(k M± , 

/i = s,i. (II. 6) 

We approximated here the scaling factors defining the 
zero-point field fluctuations with frequency-independent 
constants £ fl . The biphoton component of the wave 
function calculated in the first-order perturbation theory 
takes the form [28]: 

|*> = ^Jdt fli(t)|vac) 

^/d 2 k s ,d^d^v, (ksX ,. s;kl± ,^) 

x a) (k s± ,uj s )a f (k^,^) |vac) (II.7) 
where the probability amplitude reads: 

r-L/2 



integral expression in Eq. (|II.8[) can be given meaning- 
ful physical interpretation. Each slice of the crystal 
contributes to a biphoton amplitude \&(k s j_, w s ; ki±, Wi). 
However, the phase of this contribution changes from 
slice to slice, thus only for certain propagation directions 
the constructive interference occurs. 

The wave function given in Eq. (|II8j) describes the 
entire field emerging from the crystal into the free space. 
However, in a typical experiment the signal and idler 
photons are coupled into SMFs. For SMFs collecting light 
in the x— z plane at angles a s and on with respect to the 
z axis, the collected spatial modes can be approximated 
by gaussians centered at transverse wave vectors k, j_ = 
xw s sma s /c and kio± = — xwi sin onjc: 



it \ w > Jt \ w v ft i v 

"Mkju-L,^) = —i= cxp I — — (k,,_L - k^oJj 



fl = S,l 



(11.10) 

Here the waists w s and Wi define the spatial extent of 
the collected modes, assumed to be constant within the 
relevant spectral bandwidth. 

The wave function $(w s ,w,;) for both photons coupled 
into SMFs is given by an overlap of the wave function 
in free space ^(k^j., w s ; kjj_, wi) with the spatial profiles 
u s (k a ±,w 8 ) and Ui(ki±,Wi) of the fiber modes: 

J d 2 k sJ _ d 2 ki± u* s (k s ± , w s )u* (kjj_ , w») * (k s j_ , w s ; ko. ,Wi). 

(11.11) 

This object will be used to calculate coincidence count 
rates and spectral properties of generated photons. For 
a pump pulse amplitude in a factorable form as that in 
Eq. l|II.2j) . it will be convenient to write 

*(w.,Wi) = 4 emp K +Wi)e(w„w 4 ). (H.12) 

Here Q(w s ,Wi) can be viewed as the effective phase 
matching function for the collected modes that includes 
the geometry of the setup and the physical properties of 
the nonlinear medium. It is explicitly given by: 



dzA p (k s± +k i _ L ,w s +w J ) e(w„w 4 ) = N [ d 2 k sl _ d 2 k J± / dzA s /(k s± +k i± ) 

~ L / 2 J J-L/2 



x exp[iAk z (k sJ _,u; s ;k i± ,uj i )z} (II. 

and the Af = £p£ s £i/(2ih). The phase mismatch 

Ak z (k s ±,u! s ;ki±,u>i) is defined using the z components 
of the wave vectors of the interacting fields: 

Ak z (k s± ,w s ;k i± ,w l ) = 

= k pz (k s j_ + kjj_ , w s + w.i ) - k sz (k s ±, w s )-k iz (ki±,Wi). 

(11.9) 

These components are determined by the frequencies 
w Sl Wi and the transverse wave vectors k s j_, k^ [27|]. The 



x u* s (k s± ,w s ) <(k l± ,^)e iAMk ^ ;k ^' )z . (11.13) 

One way to simplify the above equation is to evaluate 
analytically the integral over length of the crystal, which 
yields: 

Q^(w s ,Wi) = Jd 2 k s± d 2 k l± u* s (k s± ,w s ) u*(k i± ,Wi) 

x A^ p (k s _L +ki±)sinc \^Ak z (k s ±,w s ;ki_L,Wi)j . 

(11.14) 
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However, the remaining integrals over k s j_ and ki_i_ are 
intractable analytically due to nontrivial form of phase 
mismatch Ak z and they must be performed by numeri- 
cal means. We will refer to this procedure as direct nu- 
merical integration and denote corresponding formulas 
with a superscript (D). The four-dimensional integra- 
tion task is computationally very demanding, and in the 
next two paragraphs we will present approximate meth- 
ods which reduce the computational effort to compute 
effective phase matching function 8(w s ,Wj). 



III. PARAXIAL APPROXIMATION 

The expression for the effective phase matching func- 
tion given in Eq. I|II.13P includes gaussian fiber mode 
functions u s (k s ±,ui s ) and it j (kjj_, Wj) that vanish very 
fast as the transverse wave vectors k s ± and ki± depart 
from the central observation directions k s o_L and kiox- 
This implies that little error is introduced when expand- 
ing the phase mismatch Ak z given in Eq. l|II.9j) up to the 
second order in deviations of the transverse wave vec- 
tors from k s o± and k;o±. After such an expansion the 
entire integrand in Eq. (|II.13|) takes a gaussian form in 
k s j_ and k;j_, provided that the spatial pump profile is 
gaussian as well. Consequently, one can perform all the 
integrals over transverse wave vectors analytically. This 
is a great simplification of the computational complex- 
ity of the problem, as we are now left only with a one- 
dimensional integral over z which needs to be performed 
numerically. We will call this method paraxial approxi- 
mation in analogy to the standard description of paraxial 
wave propagation in classical optics. 

It will be convenient to introduce the following nota- 
tion for the expansion of the wave vector mismatch: 

Ak z (k s± ,uj s ;k i± ,uj l ) w 
Do(w s ,Wi) + (u s ,uji)K + k t D 2 (w s ,w 1 )k, (III.l) 



where 



k = (k s _L - k s0 ±,k.i_L - k. l0± ) 7 



(111.2) 



is a four-element vector of deviations from the central 
observation directions. The Taylor expansion coefficients 
can be grouped into a scalar in the zeroth order 

D (w s ,Wj) = Ak z (k s0 ±,uj s ;k l0 ±,uj l ) (111.3) 

a vector in the first order 



i~v / s ( d s (u s ,u)i) \ 
Di w s ,Wj = , > { 



(III A) 



and a matrix in the second order: 



D2(Ws ' Wi) -' d si (co s ^), d«K, Wi ) ) (IIU) 



We wrote the vector Di and the matrix T>2 in a block 
form with entries given by: 



d M (w s , u>i) 



fdAk z dAk z 



\ 9k^ x dk^y 



(111.6) 



= ksO-L 



and 



A^ v (uj s ,uji) = 
/ d 2 Ak z 



d 2 Ak z 



d 2 Ak z 
d 2 Ak z 



(III.7) 



k s j 
k,j 



k sC 
k i0 



where /x, v = s, i. 

In order to write a compact formula for the effec- 
tive phase matching function in the paraxial approxi- 
mation, it will be helpful to represent the product of 
the fiber mode functions u*(k s j_, ui s )u*(ki± 1 u>i) and the 
pump beam profile A'^ > (k s ± + k;j_) as an exponent of a 
quadratic expression: 



u*(k s± ,Lu s )u*(k i± ,Lu i )A s I P(k s± +k i± ) = 

= exp (— Bq — B^k — k t B2«) 



(III. 



where k is a four-element vector of deviations from cen- 
tral observation directions defined in Eq. l|III.2p . The 
coefficients of the quadratic expression are a scalar 



Bo — (k S Q±_ + k M ±y 



a four-component vector 
and a 4 x 4 matrix 



p \ k s0± + k 



H0± 
iO_L 



B 2 = I ( K 



(w 2 p + w 2 )I 



(111.9) 



(111.10) 



(III.ll) 



where I denotes a two dimensional identity matrix. This 
notation allows us to write the result of four-dimensional 
gaussian integration of Eq. 13|1 over the transverse 
wave vectors as: 

n(P) , , f L/2 , AfwsWiWp 
( >(u s ,ijj l ) = I dz 



l/2 yJndctM 2 {z) 
x exp ^Mo(z) - ^Mj (z)M- 1 (z)M 1 (z)^ (111.12) 

where the superscript (P) stands for the paraxial approx- 
imation, we introduced 



Mj(z) = B 3 - izT>j 



.7 = 0,1,2 



(111.13) 
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and for notational simplicity we suppressed dependence 
on frequencies u s and Wj. The integral over the crystal 
length in Eq. (|III.12|) needs to be calculated numerically, 
which is substantially faster than direct numerical inte- 
gration of Eq. l|II.14p . It is worthwhile to note that in 
Eq. (|III.12|I the effects of spectral dispersion are fully 
taken into account, as no expansion in the signal and 
idler frequencies u> s and has been applied. As we will 
see in Sec. IVI} this makes numerical results based on the 
paraxial approximation very precise. 



IV. COSINE-GAUSSIAN APPROXIMATION 

The numerical effort to calculate the effective phase 
matching function can be reduced further at the cost of 
the accuracy. The basic idea is to replace the sine term 
appearing in Eq. 14[) by an analytically tractable ex- 
pression. Previous works [Hi EH introduced the gaussian 
approximation (GA), which approximated the sine term 
by a gaussian function, thus enabling analytical integra- 
tion. We will consider here a more general expression of 
the form: 

sine a; f=a exp(— £x 2 ) cos (C,x) = | exp(— £x 2 + i£x) + c.c. 

(IV.l) 

As seen in Fig. [21 using the parameters £ = ^ and C = \ 
yields a more accurate approximation to the sine function 
than the GA corresponding to the choice of parameters 
i = i and C = 0. 



tion directions: 




FIG. 2: (Color online) A comparison of the sines func- 
tion (solid blue line) with the cosine-gaussian approximation 
sinca; ~ exp (— x 2 /20) cos(a:/2) (circles) and the gaussian ap- 
proximation sincx ~ exp(— x 2 /5) (dashed line). 



The above observation leads us to the idea of 
cosine-gaussian approximation (CGA). Specifically, in 
Eq. (|II.14[) we replace the sine function with Eq. IjlV.lj) 
and expand the phase mismatch Afc z up to the linear 
term in transverse wave vectors around central observa- 



sinc — Afc z (k s j_,w s ;k i _ L ,tj l ) 



■ cxp 



DJk) 2 L 2 



:C(D 



k)L 



+ C.C. 



(IV.2) 



We used here parametrization introduced in Eqs. (|III.2j) - 
(|III.4|) . After inserting the above expression into 
Eq. I|II.14|) . the integrals over transverse wave vectors 
can be evaluated analytically as long as the pump and 
fiber modes are gaussian. This yields an expression for 
the effective phase matching function of the form: 

Q^{u a ,Ui) =r(u> s ,LJi)e-K u >»^cos{g(u s ,tJ i )}. (IV.3) 

The three functions appearing in the above formula are 
given by: 



2 N 



VdetK 



1 



1 



B + ]hl 2 t> 2 + ^n t k- x n 



where we defined: 



K 



N 



B x 



(£ZD + C)Di. 



(IV.4) 

(IV.6) 

(IV.7) 
(IV.8) 



For the sake of brevity we have omitted the frequency 
dependence. The expression for the effective phase 
matching function in the gaussian approximation is easily 
obtained by inserting £ = 5 and C = 0. 

In order to analyze the applicability of CGA, it is 
convenient to view the biphoton wave function given in 
Eq. (jlLllj) as an integral over k s j_ and kij_ of a prod- 
uct of two factors. The first one is the phase match- 
ing term sinc[Afc z (k s j_, ui s ; kjj_, u>i)L/2], while the sec- 
ond one, which we will call here the beam term, is a 
triple product of the pump pulse spatio-temporal pro- 
file A p (k s ± + ki±,u> s + uji) and the fiber mode profiles 
u s (k a ±,u) 8 ) and Ui(k,x,Wj). The beam term defines 
the range of transverse wave vectors and frequencies for 
which the cosine-gaussian approximation of the phase 
matching term should be accurate. This condition is sat- 
isfied when the sine argument Ak z (k s ± 7 uj s - : 'k i ±,uj i )L/2 
does not exceed approximately 3ir/2. 

Let us analyze this condition more closely. For the 
profiles assumed throughout this paper, the beam term 
takes a gaussian form: 



exp 



2c 2 



(IV.9) 



where v u 



ci>o are detunings from the central fre- takes the form [29]: 



quency and we assumed that the photons are collected 
at identical angles a s = on = a. In the exponent, we 
neglected the cross-term correlating wave vectors with 
frequencies. 

The characteristic width of the Gaussian function de- 
fines the relevant range of parameters. Thus the sum 
of the detunings is restricted by \u s + h\ < t~ x . Sim- 
ilarly the range of relevant transverse wave vectors can 
be crudely characterized by the smallest eigenvalue of the 
matrix B2 , which is equal to w s in case of symmetric cou- 
pling w s = Wi. This can be written as \k\ < w s . In the 
case of perfect phase matching for the central wave vec- 
tors k s o_L , kioj_ at the frequency ujq of the down-converted 
photons, we estimate the argument of the sine function 
expanding the wave vector mismatch Ak z up to the first 
order: 



Ak z w Di (u> , oj )k + (3(v s 



(IV.10) 



=2ujq 



dk sz I 

duj I u)=(jJq ' 



Thus we see that 



where /3 = ^ 

the CGA will be valid, if |k| and \v a + v%\ within ranges 
defined by the beam term yield the argument of the sine 
function < 'in/2. This gives: 



and 



w. 



>i|Di| 



(IV.ll) 



(IV.12) 



As the right hand sides in the above formulas are esti- 
mates, we rounded up numerical factors to simpler forms. 



V. FIGURES OF MERIT 

We will employ the computational methods presented 
in the preceding sections to analyze two parameters char- 
acterizing the usefulness of SPDC sources. The first 
one is the brightness, proportional to the probability of 
producing a fiber-coupled photon pair by a single pump 
pulse: 



R c 



(V.l) 



We will set the brightness unit by putting the multiplica- 
tive factor appearing in Eq. (|II.8|) to be \J\f\ = 1. 

The second important property of photon pairs is 
their suitability for multiphoton interference experi- 
ments. When interfering photons from independent 
sources, their spectral amplitudes cannot carry any dis- 
tinguishing information about the origin of the photons. 
This means that the biphoton wave function for each pair 
should be factorable. The degree of factorability can be 
quantified with the help of the Schmidt decomposition, 
which for the normalized wave function ^f(u! s , uji)/^/R~ c 



>R, 



= v^MKM- (V.2) 



n=0 



In the above expression, 4>^{uj s ) and <t> % n {uJi) are two or- 
thonormal sets of mode functions for the signal and the 
idler photons. The nonnegative parameters <^„ charac- 
terize the contribution of each pair of modes to the su- 
perposition. They satisfy the normalization constraint 
YlnLo = 1 and it is convenient to put them in the de- 
creasing order. Perfect factorability thus corresponds to 
the condition Co = 1- 

The degree of factorability can be quantified by the 
visibility of two-photon interference. Suppose that two 
heralded signal photons produced by identical sources 
are superposed on a 50:50 beamsplitter and the depth 
of the Hong-Ou-Mandel dip [3(J is measured. The depth 
is given by a nonnegative expression 



V 



(V.3) 



which will be called the purity parameter of a photon 
pair. In general V < 1 and the equality sign holds only 
for a factorable biphoton wave function. The purity pa- 
rameter is the inverse of cooperativity parameter intro- 
duced in Ref. [HI]. 

Typically, photon pairs are spectrally filtered in or- 
der to improve their characteristics and to lower the 
background count rates. The effects of spectral filter- 
ing can be taken into account by multiplying the two- 
photon wave function by spectral amplitude transmis- 
sions A Ai (w M ) characterizing the filters: 



(V.4) 



Note that the above substitution correctly takes into ac- 
count the decrease in count rates resulting from spectral 
filtering. We will model spectral filters using gaussian 
profiles with respective widths a s and Cj, assuming per- 
fect transmission at the peak frequency ujq: 



A^(w) = exp - 



S, I 



(V.5) 



It is worthwhile to stress that the spatial filtering imposed 
by SMFs and spectral filtering implemented with inter- 
ference filters are of different nature. The SMFs perform 
coherent filtering at the field level, i.e. add field ampli- 
tudes, while spectral filters transmit independently each 
frequency component. 

Before discussing characteristics of realistic sources, 
it is insightful to consider the limit of perfect 
phase matching, based on an assumption that 
Afc z (k s j_, uj s ] k»Xj Ui)L/2 ps over the relevant range of 
frequencies and wave vectors. This approximation means 
that we can put D = T>i = D 2 = 0, which makes the 
integrand in Eq. (|III.12j) independent of z and leads to a 
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very simple formula for the fiber-coupled biphoton wave 
function: 



W t WpW s 



■ cxp 



n 2 o (uj )w' 2 
2c 2 



(u; s a: s - LOiat)' 



x exp 



-(w s + o; i -2wo) 2 (V.6) 



where by the superscript (0) we indicated the approxima- 
tion of perfect phase matching. We also took the refrac- 
tive indices at the central frequency n (u> s ) « n (uji) w 
n o (u! ) and denoted 



1 



1 



-1/2 



(V.7) 



Let us note that the assumption Ak z L/2 « implies 
a specific geometry of the setup. First, it means that 
the pump, signal and idler beams maintain good spatial 
overlap through the entire length of the crystal. Secondly, 
the length L of the crystal must be much shorter than 
the characteristic Rayleigh range of the beams. 

The wave function given in Eq. (|V.6j) is gaussian, which 
leads to closed analytical formulas for parameters of in- 
terest. The brightness can be easily calculated to be 
equal to: 



167r 3 / 2 cL 2 



n (uj )(a s + Oi) w 2 w 2 w^ 



(V.8) 



It is instructive to analyze the scaling of the pair produc- 
tion rate in the parameters involved. The quadratic de- 
pendence on the crystal length L is a result of a coherent 
summation of the probability amplitudes of generating a 
photon pair over the entire range of —L/2 < z < L/2. 
Assuming that the waists of the pump, signal, and idler 
beams are of the same order characterized by w, the pair 
production rate scales as 1/w 3 . This scaling can be in- 
terpreted as a result of an interplay of two effects. The 
first one is the dependence of the nonlinear process on 
the transverse spatial dimension of the interacting modes. 
Suppose that the modes are confined to a transverse area 
of the order of w 2 . Then their normalization includes a 
factor 1/w for each of the modes. As the probability 
amplitude for pair generation involves an integral of a 
product of three mode functions over an area of size w 2 , 
this gives its scaling as 1/w. Squaring this result gives 
the probability of pair generation scaling as 1/w 2 . The 
second effect is the broadening of the spectrum of the 
produced photons with decreasing waists seen in the first 
exponent in Eq. (|V.6|i . which yields an additional factor 
of 1/w. 

The expression calculated in Eq. (|V.8|I enables us to 
optimize the pair production rate with respect to some 
parameters of the setup. For example, suppose that the 
waists w s and Wi of the collection modes are fixed. An 



easy calculation shows that the maximum production 
rate is achieved for the pump beam waist w p given by: 



(V.9) 



which reduces to w p = w s /2 for equal waists of collection 
modes. We will use this coupling strategy through the 
rest of the article. Note that in the case of a monochro- 
matic pump, in crude approximation of perfect phase 
mismatch the condition for optimal brightness for short 
crystal lengths takes the form w p = w s /\/2 [22l |. 

As noted in Refs. [ill . |25| . in the approximation of 
perfect phase matching the condition for spectral decor- 
relation within a photon pair is achieved when 



w p a s ai 



(V.10) 



A more general analytical condition can be derived using 
the gaussian approximation (l0| . Within this model the 
biphoton wave function takes following form: 




•(w s +Ui-2u ) /2 



(V.ll) 

Taking T(u> s ,uji) m T(ui , ui ) and expanding f(u s ,uji) up 
to the second order in frequencies around luo yields a 
gaussian expression in detunings. Spectral decorrelation 
corresponds to the vanishing cross-term (w s — u>o)(^i— u>q) 
in the exponent, which gives: 



2 r, 9 f(Us,Ui) 
T„ = Z 



duj s duji 



(V.12) 



More accurate models of the effective phase matching 
function in Eqs. (jlll. 12j) and (|IV.3j) do not yield a decor- 
relation condition in a closed analytical form. 



VI. COMPARISON 

Let us now compare computational methods intro- 
duced in the preceding sections for typical experimental 
settings. In Fig. [3] we depict the effective phase match- 
ing function <d(uj s ,uJi) for two exemplary lengths of the 
nonlinear medium calculated using direct numerical inte- 
gration, the paraxial approximation, the cosine-gaussian 
approximation and the gaussian approximation. Calcu- 
lations were carried out for a beta-barium borate crystal 
with its optical axis lying in the plane of the collected 
modes and cut at 6 C = 30° with respect to z axis. This 
corresponds to the symmetric cone half-opening angle 
equal to a = 2.2° for frequency-degenerate photons at 
780 nm. The beam waists were set to rather low values 
w s = Wi = 2w p = 70 /im to test the applicability limits 
of the paraxial approximation. 

As seen in Fig. [3j the main qualitative difference be- 
tween the computational methods is the reproduction of 
the side lobes. The impact of the side lobes on observable 
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quantities depends on the spectral width of the pump 
pulse. If the spectral bandwidth is narrower than the 
width of the central peak, then all the models can be ex- 
pected to yield similar results. Because the characteristic 
width of Q(u) a , u>i) along the axis uj s = u>i decreases with 
a longer crystal length, this regime corresponds to suf- 
ficiently narrow spectral bandwidths and short crystals. 
When leaving this regime, CGA can be expected to yield 
more accurate results in the intermediate regime com- 
pared to GA, as it reproduces correctly the lobes closest 
to the central peak. 

These predictions are confirmed by the calculation of 
the brightness R c as a function of the crystal length using 
different models, with the results shown in Fig. HJ The 
full width at half maximum of the gaussian pump pulse 
was taken equal to r^ WHM = r p Vln2 = 100 fs. The 
brightness has been calculated through two-dimensional 



numerical integration of |$(w s 



over the signal and 



the idler frequencies on a 32 x 32 square grid centered 
at u>o for the relevant frequency range where wave func- 
tion is nonzero. We have found that the further in- 
crease of grid density to 64 x 64 did not change the 
results noticeably. In the paraxial approximation, the 
effective phase matching function 6^ was evaluated at 
each point of the grid using Gauss-Kronrod quadrature 
with three-digit precision. Results based on numerical 
integration of |$(w s , uOi)\ 2 involving CGA and GA ex- 
pressions for the effective phase matching function have 
been labeled respectively as numerical CGA and numer- 
ical GA. In addition, we present results of applying a fur- 
ther simplification to CGA and GA, labeled as analytical 
CGA and analytical GA. The simplification consists in 
expanding the functions f(uj s ,uji) and g(ui s ,uji) that ap- 
pear in Eq. (|IV.3j) around the central frequency ojq up to 
the second order and replacing T(u a ,Ui) by its value at 
lu s = LOi = loq. After this expansion the squared absolute 
value of the biphoton wave function becomes a sum of 
three gaussian components and the integration over the 
frequencies ui s and u>i can be carried out analytically. 

Fig. |4] shows that for short crystals all the models give 
similar results. Furthermore, in this regime the bright- 
ness R c exhibits quadratic dependence on the crystal 
length, which agrees with Eq. (|V.8|) derived under the as- 
sumption of perfect phase matching. As expected, with 
an increasing crystal length the GA model departs earlier 
from the numerical results than the CGA model. 

A more thorough way to compare the paraxial approx- 
imation with direct numerical integration is to evaluate 
two quantities: the scalar product between the normal- 
ized biphoton wave functions $' p ^ and <i>^ D ^ obtained 
using both methods and the ratio of the corresponding 
pair production rates Rc /Rc - We carried out these 
calculations in an unfavorable regime of a long crystal 
L = 2 mm, ultrashort pump pulses r PWHM = 20 fs, and 



strong focusing w s 



2u>„ 



40/jm. We found 



that both the quantities differed from one by less than 
10~ 3 . It should be noted that the computational effort 
required by the paraxial approximation was reduced in 
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FIG. 3: The effective phase matching function Q(u) s ,uJi) cal- 
culated using (a, b) direct numerical integration; (c, d) parax- 
ial approximation; (e, f) cosine-gaussian approximation and 
(g, h) gaussian approximation for the crystal length (a, c, e, 
g) L = 100 ^m and (b, d, f, h) L = 1 mm. The pump and 
collecting beam waists were set to w 3 = Wi = 2w p = 70 fim. 
The angular frequencies io s and are labelled with the cor- 
responding wavelengths. 
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FIG. 4: (Color online) The source brightness R c calculated 
using different numerical methods, specified in the inset, as 
a function of the crystal length L, for the waists w 3 = Wi = 



2w p = 70 (Jm and the pump pulse duration r p 



100 fs. 



our calculations by 
ical integration. 



10 4 compared to the direct numer- 



Finally, let us analyze the coincidence count rate R c as 
a function of the pump beam waist w p and the fiber mode 
waists in a symmetric setup, when w s = Wi. In Fig.[5]we 
depict results obtained using the paraxial approximation 
for two exemplary lengths of the crystal. It is seen that 
for a fixed waist of the fiber modes the brightness has a 
well pronounced maximum in w p . This maximum is lo- 
cated to a good approximation at w p = w s /2, which is in 
an agreement with the result derived within the elemen- 
tary model of perfect phase matching in Eq. I|V.9|1 . This 
motivated the choice of w s — u>i — 2w p in the presented 
examples. 
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FIG. 5: (Color online) The natural logarithm of the brightness 
In R a as a function of pump beam waist w p and fiber mode 
waists w s = Wi for the crystal length (a) L = 1 mm and 
(b) L — 100 /im. The dashed (red online) lines depict the 
condition w p = w s /2 specified in Eq. (|V.9[) . 



VII. SPECTRALLY UNCORRELATED PAIRS 

A necessary condition for high-visibility multipho- 
ton interference is the lack of distinguishing informa- 
tion about the origin of the photons, which means that 
each the photon should be prepared in an identical pure 
wavepacket. The most obvious way to achieve this 
regime is to insert interference filters whose bandwidth 
is smaller than the characteristic scale of spectral cor- 
relations within photon pairs. An intriguing alternative 
has been presented in Ref. which proposed to re- 
move spectral correlations by exploiting geometric effects 
in SPDC. The purity of the produced photons needs to 
be analyzed in conjunction with other characteristics of 
the source, such as the pair production rate. In this sec- 
tion we will employ our computation tools to compare 
properties of spectrally decorrelated pairs generated by 
different methods. 

Let us first analyze the geometric approach of Ref. [ll[. 
The underlying physics can be understood intuitively 
by looking at the biphoton wave function in the perfect 
phase matching approximation given by Eq. (|V,6|) . The 
spectral pump amplitude introduces anticorrelations be- 
tween frequencies of the down-converted photons, while 
the pump beam waist and emission angles define the de- 
gree of positive correlations. By balancing these two ef- 
fects one can obtain a factorable biphoton wave func- 
tion. More generally, without the approximation of per- 
fect phase matching, one needs to analyze correlations in- 
troduced by the function 0(w s ,o;,) defined in Eq. 12j) 
combined with the spectral pump amplitude. As the non- 
linear medium we considered a BBO crystal in the same 
configuration as discussed in Sec. EH As the basic tool, 
we chose the paraxial method developed in Scc. lIIII duc to 
its high precision and computational effectiveness. In or- 
der to evaluate the purity parameter V measuring degree 
of spectral correlations, the approach presented by Law 
et al. [29] was used. The method is based on the sin- 
gular value decomposition of the matrix representation 
of the biphoton wave function $(w s ,Wj) on a sufficiently 
fine discrete grid. The normalized singular values are ap- 
proximations of Schmidt coefficients <; n and as such are 
used to evaluate purity parameter V . We found it suf- 
ficient to take the grid 32 x 32. Further increase of the 
grid density did not make any noticeable difference. 

In Fig. [6](a,b) we present the purity parameter for two 
typical lengths of the crystal as a function of the pulse 
duration T p WKM and the collecting mode waist w s . We 
assumed that the waists of the fiber modes and pump 
beam are w s — u>i — 2w p , which is motivated by the 
results presented in Fig. \E[ The contour plots exhibit 
a clear relation between r p WHM and w s that leads to 
minimized spectral correlations between photons. For 
a comparison, Fig. [U[a) and (b) depict also the purity 
condition derived in Eq. l|V.12p using the GA model, 
as well as the predictions of the perfect phase match- 
ing approximation given in Eq. (|V.1Q|1 . It is seen that 
for the shorter crystal length L — 100 ^m the simple 
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(a) 




FIG. 6: Contour plots of (a, b) the purity parameter V and 
(c, d) the brightness R c as a function of the pump pulse du- 



ration T„ 



and the collected mode waists w s = Wi. The 



crystal thickness is (a, c) L = 1 mm and (b,d) L — 0.1 mm. 
The solid and dashed lines in (a, b) correspond to factorabil- 
ity conditions given respectively in Eqs. (|"V. 12[) and (|V.10I) . 
The grey areas in (c, d) mark the regions where the purity 
parameter is V > 0.99. 



analytical formula of Eq. IjV.lOp gives accurate results. 
This is because the spectral anticorrelations are predom- 
inantly defined by the pump bandwidth rather than the 
phase matching of the crystal. This is no longer valid 
for the length L = 1 mm, where the effective bandwidth 
of the down-conversion process becomes strongly affected 
by the phase matching. These observations are consistent 
with results presented in Fig. |4j for L = 100 /im the pair 
production rate is accurately given by the perfect phase 
matching approximation, while for L = 1 mm effects of 
finite phase matching bandwidth are clearly seen. 

The relation between the collecting mode waist w s = 
Wi and the pump pulse duration rJ WHM that leads to 
minimized spectral correlations gives us some flexibility 
to optimize the source with respect to other parame- 
ters. In Fig. [UJc,d) we present the source brightness R c 
as a function of w s and r p . Note that in our calcula- 
tions we constrain the pump beam waist by imposing 
w s = Wi = 2w p . It is seen that R c can be increased by 
reducing the fiber mode waist w s . As Figs. (He) andEJd) 
depict the pair production rate in the same units, we can 
compare the brightness for the two crystal lengths. As- 
suming that we have no restrictions on the pump pulse 
duration, a shorter crystal can produce more uncorre- 
cted photon pairs. This is because for L = 1 mm 
stronger spectral anticorrelations overwhelm the benefit 



of a longer nonlinear medium. However, in a realistic sit- 
uation there is usually a technical minimum on the pump 
pulse duration. For concreteness, let us assume it to be 
t fwhm = 100 fe An inspection f Fig. [g] snows that 

under the condition of nearly ideal decorrelation defined 
by the value of the purity parameter V ~ 0.99 higher 
brightness, approximately equal to R c « 0.046, is ob- 
tained when the fiber mode waist is w s ~ 1 mm and the 
crystal length L = 1 mm. We found that for even longer 
crystals decorrelation can be reached only using longer, 
less focused pump pulses, which lowers the source bright- 
ness. 

These limitations raise the question whether a more 
efficient strategy may rely on collecting tightly focused 
modes and removing spectral correlations with interfer- 
ence filters. Let us consider the same pump pulse dura- 
tion Tp WHM = 100 fs and crystal length L — 1 mm as 
before, but tighten the fiber mode waists to w s = 100 ^m. 
The result is significantly increased brightness, but at the 
cost of introducing spectral correlations. The effects of 
inserting interference filters into such a setup are shown 
in Fig. [7l where we depict the brightness R c and the 
purity parameter V as a function of the spectral filter 
bandwidth. It is seen that for the bandwidth a ~ 2.6 nm 
the purity parameter reaches the value V « 0.99, while 
the brightness is R c ~ 3.8, which is significantly higher 
than before. Thus the benefit of increased brightness is 
retained despite spectral filtering. 
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FIG. 7: (Color online) The brightness R c (dashed blue line, 
left vertical scale) and the purity parameter V (solid red line, 
right vertical scale) as a function of the spectral filter band- 
width a — a a = Oi for a crystal length L = 1mm, beam 
2w p = 100/^m and the pump pulse dura- 



waists w £ 
tion r v = 



= Wi 

100 fs. 



In order to gain more insight into the trade-off between 
the source brightness and spectral correlations, we calcu- 
lated the maximum filter bandwidth that gives the purity 
V ~ 0.99 for a range of pump beam waists w p , while keep- 
ing other parameters of the setup identical as in previous 
examples. The results are shown in Fig.[8l It is seen that 
the filter bandwidth across the analyzed range does not 
deviate significantly from the value a = 2.7 nm, while 
the brightness increases substantially with tighter focus- 
ing. This can be explained by the fact that the spectral 
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filter bandwidth is defined by the requirement to remove 
frequency anticorrelations which depend primarily on the 
crystal length and the pump pulse duration rather than 
the beam waist. 
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FIG. 8: (Color online) The brightness R c (solid red line, left 
vertical scale) as a function of the collected mode waist w s 
obtained for the maximum filter bandwidth (dashed blue line, 
right vertical scale) which yields the purity parameter above 
V > 0.99. Other setup parameters are identical as in Fig. [7] 



VIII. CONCLUSIONS 

In this paper we introduced and utilized approximate 
methods that alleviate the numerical load necessary to 
model SPDC sources while retaining the accuracy of the 
results in physically relevant regimes. Our approach 
was based on an observation that optical fibers collect- 
ing photons effectively define a relatively narrow range 
of wave vectors that needs to be included in calcula- 
tions. This justified applying the paraxial approxima- 
tion, which made a substantial portion of the problem 
tractable analytically and significantly reduced the re- 
maining numerical effort. The paraxial approximation 
can be also combined with a simplification of the two- 
photon wave function to an analytically manageable form 
that led to closed formulas. We exploited these strate- 
gies to analyze performance parameters that characterize 
the usefulness of SPDC sources for quantum information 
applications, such as the pair production rate and the 
spectral purity parameter that is critical in multiphoton 
interference experiments involving multiple sources. 

The choice of a computation method depends on the 
range of the setup parameters. The most difficult regime 
to deal with is that of very broadband, tightly focused 
pump pulses and long crystals. It is then necessary to 
include with high precision the phase matching function 
over a wide range of frequencies and transverse wave 
vectors. The most universal method is then direct nu- 
merical integration, which however requires tremendous 
computational effort. In practical situations, the parax- 
ial approximation, delivers highly accurate results with 



significantly reduced numerical load for typical setup pa- 
rameters. The validity of the paraxial approximation 
can be checked with a relatively low effort by compar- 
ing it with direct numerical integration only at the edges 
of the region of interest that correspond to most unfa- 
vorable cases. Such a confirmation allows one to apply 
the paraxial approximation throughout the entire region 
of interest reducing the overall computational cost. For 
examples studied in Sec. IVII[ the paraxial approxima- 
tion has been verified to yield results that did not differ 
by more than few percent from direct numerical inte- 
gration. In more restricted scenarios, one may consider 
using the cosine-gaussian approximation, which extends 
the validity of the previously used gaussian approxima- 
tion. Results obtained with these methods can be used 
as a starting point for designing source characteristics 
with more elaborate and precise tools. We also discussed 
a crude approximation of perfect phase matching, which 
gives simple, qualitative insights into the roles played by 
various source parameters. 

The numerical methods presented in this work can 
be used to analyze various aspects of down-conversion 
sources that are relevant to experimental implementa- 
tions of quantum information processing protocols. We 
discussed here spectral decorrelation, which is a necessary 
condition to achieve high- visibility multiphoton interfer- 
ence between independent sources, in connection with 
the pair production rate. For exemplary settings chosen 
for the analysis, we found that spectral filtering com- 
bined with tight focusing of the pump beam can deliver 
higher brightness than balancing the spectral correlations 
using the geometry of the setup. The paraxial approx- 
imation can be also extended to analyze properties of 
an individual photon generated in the down-conversion 
process, with traced out degrees of freedom of the con- 
jugate photon. This approach has been successfully ap- 
plied to model the results of a measurement of the single- 
photon density matrix in the spectral domain reported 
in Ref. (24|. Theoretical details of this work will be pre- 
sented elsewhere Furthermore, the single photon 
count rates allow us to calculate the heralding efficiency, 
defined as the ratio of the pair production rate to the 
count rate on the trigger detector. This is another im- 
portant parameter characterizing the usefulness of down- 
conversion sources [T3|, that can be efficiently submitted 
to numerical optimization using paraxial approximation. 
We aim to make this a subject of a separate publication. 
The numerical results presented in this paper have been 
obtained using a Mathematica code which can be down- 
loaded from [33|. 
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